Quantum Interference in Superconducting Wire Networks and Josephson Junction 
Arrays: Analytical Approach based on Multiple-Loop Aharonov-Bohm Feynman 

Path-Integrals 
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We investigate analytically and numerically the mean-field superconducting-normal phase bound- 
aries of two-dimensional superconducting wire networks and Josephson junction arrays immersed in 
a transverse magnetic field. The geometries we consider include square, honeycomb, triangular, and 
kagome lattices. Our approach is based on an analytical study of multiple-loop Aharonov-Bohm 
effects: the quantum interference between different electron closed paths where each one of them en- 
closes a net magnetic flux. Specifically, we compute exactly the sums of magnetic phase factors, i.e., 
the lattice path integrals, on all closed lattice paths of different lengths. A very large number, e.g., 
up to 10 81 for the square lattice, exact lattice path integrals are obtained. Analytic results of these 
lattice path integrals then enable us to obtain the resistive transition temperature as a continuous 
function of the field. In particular, we can analyze measurable effects on the superconducting transi- 
tion temperature, T C (B), as a function of the magnetic filed B, originating from electron trajectories 
over loops of various lengths. In addition to systematically deriving previously observed features, 
and understanding the physical origin of the dips in T C (B) as a result of multiple-loop quantum 
interference effects, we also find novel results. In particular, we explicitly derive the self-similarity 
in the phase diagram of square networks. Our approach allows us to analyze the complex structure 
present in the phase boundaries from the viewpoint of quantum interference effects due to the elec- 
tron motion on the underlying lattices. The physical origin of the structures in the phase diagrams 
is derived in terms of the size of regions of the lattice explored by the electrons. Namely, the larger 
the region of the sample the electrons can explore (and thus the larger the number of paths the 
electron can take), the finer and sharper structure appears in the phase boundary. Our results for 
kagome and honeycomb lattices compare very well with recent experimental measurements by Ciao, 
Hues, Chaikin, Higgins, Bhattacharya, and Spencer (Phys. Rev. B, companion paper). 

PACS numbers: 74.50. +r 



I. INTRODUCTION 



A. Physics of the Phase Diagram 
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When immersed in an externally applied magnetic 
field, superconducting networks^ made of thin wires, 
proximity-effect junctions, and tunnel junctions ex- 
hibit complex and interesting forms of phase dia- 
grams. These superconducting networks have been stud-, 
ied in various, kinds of geometries, including .simplcEl 
and complexoa periodic lattices, regular Jractalsjj bond- 
percolation networks J3 disordered arraysu and quasiperi- 
odic lattices.Lnlil The rich structure present in the 
resistive transition temperature as a function of the 
magnetic field, namely, the superconducting-normal 
phase diagram, has a rich structure that has been 
the subject c£_j/arious experimental and theoretical 
investigations .Eni3Ti3 



The rich structure in the phase diagram is essentially a 
result of quantum interference effect or frustration due to 
the magnetic field and the built-in multi-connectedness 
of the networks. The magnetic fluxes through the cells of 
various areas, measured in units of the superconducting 
flux quantum $o = hc/2e, are useful parameters to char- 
acterize the interference effect. At zero magnetic field, 
the quantum interference effect is absent, and therefore 
the resistive transition temperature should have a peak. 
Also, due to gauge invariance, physical quantities should 
be periodic functions of the cell fluxes, with a period of 
$o- These arguments qualitatively explain the apparent 
periodic or quasi-periodic structures observed in phase 
diagrams of networks of various geometries. 

To gain a quantitative description of the phase dia- 
grams, we employ the mean-field theory which is very ef- 
fective in serving such a purpose. For wire networks, the 
mean-field expression is given by the Landau-Ginsburg 
equation expressed in terms of the order parameters at 



1 



the nodes Jiia For a jiiaetion array, one has a set of self- 



consistent equationsti-ll__l for the thermally averaged pair 
wavefunctions of the grains. Such equations are lin- 
earized near the transition point, and the highest tem- 
perature at which a non-trivial solution first appears is 
identified as the transition temperature. Therefore, one 
is left to find the top spectral edge of eigenvalue prob- 
lems. The equations for a junction array can be mapped 
into a tight-binding Schrodinger problem for an electron 
hopping on a lattice immersed in a magnetic field. The 
equations for a wire network are in general more difficult 
to solve, because the eigenvalue appears in a non-linear 
way. 

Numerical results have been obtained for phase dia- 
grams of networks of various geometries. All of them 
compare very well with the corresponding experimental 
data; the locations of the peaks of various sizes are cor- 
rectly predicted and the relative heights of the peaks are 
also reproduced with occasional small deviations. The 
success of mean-field theoryliSlLj suggests that much of 
the frustration effect in a statistical problem can be ac- 
counted for in terms of quantum interference effect of 
linear wave mechanics. 



B. Many-loops generalization of the standard 
Aharonov-Bohm effect 

In this paper, we systematically investigate the field- 
dependent superconducting-normal phase for a variety of 
two-dimensional superconducting networks. The basis of 
our approach is the analytic study of electron quantum 
interference effects originating from sums over magnetic 
phase factors on closed lattice paths. The sums of these 
phase factors, called lattice path integrals, are many- 
loop generalizations of the standard one-loop Aharonov- 
Bohm- type argument, where the electron wave function 
picks up a phase factor e 1 * each time it goes around a 
closed loop enclosing a net flux $. 

We compute analytically the lattice path integrals up 
to very long lengths for various types of lattices. These 
lattice path integrals contain the quantum interference of 
enormous numbers of closed paths. Through an iterative 
approach, these results then-enable us to obtain the corre- 
sponding phase boundariesEjtJ as continuous functions 
of the strength of the applied field. This method pro- 
vides a systematic approximation scheme, through finite 
truncations, for the spectral edges of eigenvalue problems 
from which our mean-field phase diagrams can be com- 
puted. Thus, we can gain considerable theoretical insight 
into the physical origin of the structure in the phase di- 
agrams. This approach also enables us to analyze the 
structure of the phase boundaries from the viewpoint of 
the geometric features of the networks. We apply this 
approach to study the phase boundaries of square, hon- 
eycomb, triangular, and kagome lattices. Our studies 
provide complete and detailed analysis of the relation- 



ship between the phase diagram structures and the cor- 
responding network geometries. 



C. Organization of the paper 

This paper is organized as follows. In Sec. II, we de- 
scribe the general formulation of our approach to the 
determination of phase diagrams for a variety of peri- 
odic superconducting networks. To illustrate our calcu- 
lational scheme, we first compute the Little-Parks oscil- 
latory phase boundary of a single superconducting loop 
in Sec. III. In Sec. IV, we apply this approach to the su- 
perconducting square network. We devote Sec. V to the 
discussion of a very important and interesting feature 
observed in the phase boundary of the square network, 
namely, the self-similarity. The superconducting honey- 
comb, triangular, and kagome networks are studied based 
on the same approach, respectively, in sections VI, VII, 
and VIII. In Sec. IX, we discuss some general trends in 
the application of this approach to these types of net- 
works studied above. Comparisons of the phase bound- 
aries between a single superconducting loop and the cor- 
responding superconducting network are also made. Fur- 
thermore, we present a brief discussion on the relation- 
ship between our approach and other related methods. 
In Sec. X, we compare the phase boundaries of honey- 
comb and kagome lattices. The last section summarizes 
our results. 



II. GENERAL FORMALISM 

The physics of T C (B), the superconducting-normal 
phase boundary as a function of the field B, is determined 
by the electronic kinetic energy because the applied fiel 
induces a diamagnetic current in the superconductor 
This current (proportional to the velocity) determines 
the kinetic energy of the system. In other words, the ki- 
netic energy can be written in terms of the temperature 
as 



2m 



2to*£(T) z 



T C (B) - T c (0), 



where, for any superconductor, m* is twice the electron 
mass, and 



y/l-T c (B)/T c (0Y 



is the temperature-dependent coherence length. The 
problem of obtaining T C (B) is then mapped to that of 
finding the spectral edges of tight-binding electrons on 
the corresponding lattice. Thus, assuming a unit hopping 
integral between adjacent sites, we consider the following 
Hamiltonian 



2 



H 



E 



c\cj exp(iAij), 



(1) 



which describes the kinetic energy of electrons hopping 
on a discrete lattice subject to a perpendicular magnetic 
field. Here (ij) refers to nearest-neighbor sites and the 
magnetic phase 



Aij = 2vr / A • dl 



is 2tt times the line integral of the vector potential, A, 
along the bond from j to i in units of the $ = hc/2e. 



A. Sums over closed paths 



The lattice path integral, is defined as 



W = E 

All closed lattice paths 7 of length I 



(2) 



By closed paths of length I we mean the paths starting 
and ending at the same site after traversing / steps on 
the lattice and <I> 7 is the sum over phases of the bonds 
on the path 7. Let denote a localized electron state 
centered at site i. It is not difficult to notice that /z; corre- 
sponds precisely to the quantum mechanical expectation 
value ( 1 &j|.ff |\E r j), which summarizes the contribution to 
the electron kinetic energy of all closed paths of Z-steps. 
The physical meaning of the lattice path integral 



thus becomes clear. The Hamiltonian H is applied I times 
to the initial state resulting in the new state 

|*/> = -ff'|*i> 

located at the end of the path traversing I lattice bonds. 
Because of the presence of a magnetic field, a magnetic 
phase factor e lAij is acquired by an electron when hop- 
ping from j to the adjacent site i. The lattice path inte- 
gral (ii is nonzero only when the path ends at the starting 
site. In other words, fil is the sum of the contributions 
from all closed paths of I steps starting and ending at the 
same site, each one weighted by its corresponding phase 
factor e 1 *^ where 



$7 

— = net flux enclosed by the closed path 7. 

Z7T 



B. Quantum Interference 

It is important to stress tha t ffi . f . depends crucially on 
the traveling route of the path.liirlia For instance, <I> 7 will 
be positive (negative) by traversing a polygon loop coun- 
terclockwise (clockwise). Therefore, quantum interfer- 
ence information contained in m arises because the phase 



factors of different closed paths, including those from all 
kinds of distinct loops and separate contributions from 
the same loop, interfere with each other. Sometimes, the 
phases corresponding to subJpops of a main path cancel. 

To analytically computeO_!rEj the lattice path integrals 
Hi is in general a difficult task since hi involves an enor- 
mous number of different paths (growing rapidly when I 
increases) , each one determined by its corresponding net 
magnetic phase factor. We have developed systematic 
and efficient methods to calculate the lattice path inte- 
grals for a number of distinct lattices. The techniques 
involve successively iterating the constructed recursion 
relations and exploiting the symmetries of the underlying 
lattices. The technical details of the implementation will 
be presented elsewhere. Below we will only list the first 
few calculated lattice path integrals in relevant places. 
Results for the lattice path integrals of larger I will not 
be presented due to their lengthy expressions, but will be 
used in some of our calculations. 

In summary, the lattice path integrals summarize the 
electron quantum interference effects originating from 
sums over magnetic phase factors on closed lattice paths. 
The sums of these phase factors, the lattice path inte- 
grals, are many-loop generalizations of the standard one- 
loop Aharonov-Bohm-type argument, where the electron 
wave function picks up a phase factor e l * each time it 
goes around a closed loop enclosing a net flux <I>. 



C. Computation of the energy eigenvalues from 
lattice path integrals 

We now outline the scheme for obtaining the eigen- 
values from the calculated lattice path integrals. Let us 
apply the Hamiltonian to the starting state 

l^i) = I*,), 

which is a localized state centered at an arbitrary site i 
on the lattice, and perform the following expansions: 



and for n > 1 



H\lp n ) = fcnlV'n-l) + 0>n\i>n) + b n+ i\i> n+1 ) . 

The Hamiltonian matrix in the basis \ip n ) is obviously in a 
real tridiagonal form. Each new state in this method ex- 
pands outward by one more step from the site where the 
starting state is located. Note that the a n 's and b n+ i's 
are gauge-invariant quantities. Through these parame- 
ters we can construct the truncated Hamiltonian matri- 
ces, H^ n \ which is the nth-order approximation to H. 
For instance, 
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b 2 a 2 
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H (3) 



and so on. The quantity we desire, i.e., the top spectral 
edge can then be obtained by solving the eigenvalues of 
and will be designated by T c , which is the nth- 
order approximant to the phase boundary. This scheme 
is useful because finite truncations give good approxima- 
tions to T C (B). 

The coefficients a„'s and 6 n +i's can be exactly ex- 
pressed in terms of the lattice path integrals in a system- 
atic manner, which will be presented below, respectively, 
for the bipartite and non-bipartite lattices. In general, 
given the lattice path integrals up to the order H2L-1, 
which contains information on the quantum interference 
effects due to closed paths of 2L — 1 steps, we can ob- 
tain the coefficients up to and b^. Thus, the Lth-order 
truncation of the Hamiltonian matrix can be constructed, 
and subsequently can be obtained. 



1. For bipartite lattices 

We first discuss the case for bipartite lattices where the 
lattice path integrals of odd number steps are identically 
zero, i.e., 



It is evident that 



M2(+i = 0. 



a n = 



for any n. To compute the b n +is, we define an auxiliary 
matrix B with the first row elements given by 

B\.i = [l 2 l. 

The other rows are evaluated by using only one immedi- 
ate predecessor row. Namely, for k > 2 and I > 1 



B. 



k,l 



Bk-id+i 
Bk-1,1 



1-1 

i=0 



(3) 



where 



B n n = 1 



for n > 1. The 6 n +i's are obtained from the elements of 
first columns of the matrix B as 



b n +i — v 



(4) 



Below we explicitly express the first few b n+ is in terms 
of the lattice path integrals noting that U2 is always equal 
to z, the coordination number of the lattice: 

b2 = VM2 = Vz, 

ki = \ — : ■ 



2/14Z + z 3 « 4 — z 2 



/U4 — -2 Z 

These expressions are applicable to any type of bipartite 
lattice. 

It is worthwhile to point out that the number of ele- 
ments on a specific row is always less than that on the 
immediate predecessor row by one. For instance, for a 
specific fc, if the matrix elements run from Bk.i to -B/c,z, 
the elements in the next row run from Bfc+1,1 to Bk+i,i-i- 
Therefore, given the lattice path integrals up to u 2 l , the 
matrix B consist of L rows. The Lth (last) row has only 
one element Bt, \ from which we can deduce 6l+i. It is 

clear now that the highest-order approximation Tc L+1 ^ to 
the phase boundary can be obtained from /j, 2 , "4> • • • , M2L- 



2. For non-bipartite lattices 

Turning to the non-bipartite lattice case, we now define 
an auxiliary matrix N with the first row elements given 

by 

The other rows are evaluated by using only one immedi- 
ate predecessor row. Namely, for k > 2 and / > 1 



1-1 



Nk,i = N ^- N ^ ; f^ £ Nk,%Nk—i,i—i, 

»=o 

(5) 



where N n fl = 1 for n > 1. The a„'s and 6 n +i's are ob- 
tained from the elements of the first and second columns 
as 



a n = N nA 



and 



bn+l = ^Wn,2 - JV£i- 



(6) 



(7) 



Below we explicitly express the first few o„'s and 6„ + i's 
in terms of the lattice path integrals: 

01 = 0, 

M3 

a 2 = — , 

z 

— H3Z — Z 4 



4 



and 



The above expressions are valid for any type of non- 
bipartite lattice. 

It is worth stressing that the number of elements on 
a specific row is always less than that on the immediate 
predecessor row by two. For instance, for a specific k if 
the matrix elements runs from Nk,i to Nk.i, the elements 
in the next row runs from Nk+1,1 to Nk+ij-2- Therefore, 
given the lattice path integrals up to [i 2 l+i, the matrix iV 
consist of L + 1 rows. The Lth row has only three element 
Nl,i, Nl^, and where can be obtained from 

Nl,2, and Nl^. The (L + l)th (last) row has only one 
element Nl+i,i from which we can deduce cll+i- It is 
clear now that the highest-order approximation Tq L+1 ^ 
to the phase boundary can be obtained from [11,(12,..., 
and [12L+1- 



III. SIMPLE ILLUSTRATION: A SINGLE 
SUPERCONDUCTING LOOP 

Before we study the lattice cases, we apply the for- 
malism described above to three simple single-cell cases. 
Namely, we calculate, respectively, the transition tem- 
perature of a single superconducting loop in the shape 
of a square, a hexagon, and a triangle. Exact solutions 
of the phase boundaries can be obtained for these simple 
cases. For all of these, $ = 0/2-7T stands for the magnetic 
flux through these elementary cells, in units of $0. 

The lattice path integrals, [n, now correspond to the 
sums over all closed paths of I steps on a single cell. 
Closed-form results for the lattice path integrals are de- 
rived. They are, respectively, 

[1/2] 

Hi? =Cf l + 2j2cfL 2k cos(kcj ) ) 



on a square, 



fc=i 



[i/3] 



^ ) ^Q 2 ' + 2^Q 2 l 3fe cos(M 



k=l 



on a hexagon, and 



[1/3] 



H$=C? l + 2j2C?U k cos(2k<P), 
fe=i 

[(/-l)/3] 

^+1=2 E C^_ lC os[(2fc + l) 



fe=0 



on a triangle. Here 



a 



ml 



is the binomial coefficient, and the notation [x] means the 
largest integer equal to or smaller than x. Through these 
results for the lattice path integrals, it is straightforward 
to compute the parameters a„'s and 6 n +i's. In fact, for 
these small simple systems, the iterative process termi- 
nates very quickly. In other words, the parameters a n 's 
and b n+ i's become identically zero after a few iterations. 
Hence, the corresponding exact tridiagonal Hamiltonian 
matrices can be readily constructed. 



Square loop 

Denoting the tridiagonal Hamiltonian matrix for the 
square loop by H Sl we find that 




1 

( ( •:■ 





1 











cos(|) 





(*) 





sin (§) 







sin (|) 






which is obtained by using only [12, [14, and [i 6 . A closed- 
form expression for the top eigenvalue of H s can be easily 
obtained 



T c {cj)) = J2 + 2cos 



Hexagonal loop 



Similarly, denoting the tridiagonal Hamiltonian matrix 
for the hexagon loop by Hh, we find that 



V2 
V2 
1 









^/l+cosO) 

^/l+cosf^) A /l-cos( 

y/i- cos {4>) 

1 



which is obtained by using only (12, [14, [i§, [is, and [i w . 
Let j be an integer, the top eigenvalue of Hh can be 
expressed as follows: 



2 + 2 cos (| + for _| +3j <^<-i+3j 

2 cos (I) for -i +6 j<A<i +6j 

2 + 2cos(|- 2 ^) for i +3 j<^:<|+3j 

-2cos(f) for | +6 j<A<i +6j 



n!(m — n)! 



•5 



Triangular loop 

Denoting the tridiagonal Hamiltonian matrix for the 
triangle loop by H t , we find that 



H t = 



V2 
\/2 cos (0) 
|sin(0)| 





|sin(0)| 

- COS {(f)) 



which is obtained by using only \i\ through /15. The top 
eigenvalue of H t can be expressed as follows: 



T c {4>) 



2 cos 1 



2tt 
3 



for -§ +3 j<^<-i+3j 

for -5+3j<^:<i+3j 

2cos(§-^l) for i+3j<^<|+3j 



2 cos ( § 



In Fig. 1, we plot the superconducting transition tem- 
perature, AT C ($) = T c (0) - T c ($) = 2 - T c ($), respec- 
tively, of a square loop, a hexagon loop, and a triangle 
loop for —2 < $ < 2. It is evident that these phase 
diagrams are qualitatively identical. Also, the AT C (<I>) 
shown are periodic functions of $ and the period of the 
oscillation in the flux is equal to $o- As expected, AT C (<1>) 
have their minima at <f> = j$o and their maxima at 
$ = j$ /2. 

It is interesting to note that AT C ($) has the largest 
magnitude for the triangular loop and the smallest for the 
hexagonal loop. It will be seen in Sec. X that this one- 
loop general behavior carries over to the network cases, 
in spite of the distinctive differences in the fine structure 
of their phase boundaries. These results are consistent 
with the ones obtained numerically in Rcf. 2. 



IV. SQUARE LATTICE 

For the square lattice, we denote the lattice path inte- 
grals by S2/ . In other words, s 2 ; is the exact sum of the 
phase factors of all 2Z-step closed paths on the square 
lattice. Below 0/27T corresponds to the magnetic flux 
through an elementary square plaquette, i.e., 

2tt 

Throughout this paper, c denotes the lattice constant of 
all the lattices considered in this work. The results for 
s 2 , S4,,..., S12 are: 

82 = 4, 

Si = 28 + 8 cos0, 
s 6 = 232 + 144 cos + 24 cos 20, 
s 8 = 2156 + 2016 cos + 616 cos 20 + 96 cos 30 
+16 cos 40, 

sio = 21944+ 26320 cos 0+ 11080 cos 20 + 3120 cos 30 
+840 cos 40 + 160 cos 50 + 40 cos 60, 



S12 = 240280 + 337560 cos 0+ 174384 cos 20 

+67256 cos 30 + 23928 cos 40 + 7272 cos 50 
+2400 cos 60 + 528 cos 70 + 144 cos 80 + 24 cos 90. 

We have computed the lattice path integrals for the 
square lattice up to S138, which are obtained by exactly 
summing up ~ 10 81 closed paths. The first few lattice 
path-integrals can be quickly obtained analytically by 
hand. We have used Maple symbolic manipulation soft- 
ware to obtain lattice path integrals of longer lengths. 
For these, it is convenient to optimize the algorithm by 
exploiting the symmetries of the problem. These calcu- 
lated lattice path integrals s 2 ; 's have enabled us to obtain 
the phase boundary up to Tc 7n \(f>). 

It is instructive to explain how the first few lattice path 
integrals are obtained. This will also clarify their physical 
meaning. Since there is no path of one step for returning 
an electron to its initial site, si is always equal to zero. 
Indeed, all lattice path integrals s 2 i+i involving an odd 
number of steps are equal to zero. Now let us compute 
the next lattice path integral, with two-steps. There are 
four closed paths of two steps each [retracing each other 
on one bond (•<-»•), where the dot • indicates the initial 
site], thus 



S2 



= 4 • =4 e m4> = 4 = 



where z is the coordination number of the lattice. 

There are 28 closed paths of four steps each: four re- 
tracing twice on one bond (■«-►); twelve starting from a 
site connecting two adjacent bonds and retracing once on 
each bond (■«-»•■«-►); and twelve moving two bonds away 

and then two bonds back to the original site (•—►—►)■ 
Since all of them enclose no area (i.e., no flux), then 



no flux 
S 4 



+ 12- 



12 



28. 



Among the 4-step closed paths, 8 of them enclose ad- 
jacent square cells (4 counterclockwise and 4 clockwise) 
contributing 



4e l0 + 4e~^ = 



<COS( 



to Si. Thus it follows that S4 = 28 + 8 cos0. Higher- 
order integrals S21 can be similarly constructed. 

It is straightforward to compute the non-zero parame- 
ters b n from the obtained results for S21 ■ The correspond- 
ing truncated Hamiltonians, H^ n \ can then be readily 
constructed. For instance, the second-order truncation 
of the Hamiltonian is 



2 
2 



Its corresponding top eigenvalue is Tc 2 \<j)) — 2, which 
does not depend on 0. This is understandable from the 
fact that the shortest length for a closed path on the 
square lattice to enclose the magnetic flux is for I = 4 



6 



while only contains elements derived from /j,2- The 
third-order truncation of the Hamiltonian is 







2 __ 

2 V3 + 2 cos ( 

V3 + 2 cos (j) 



Its corresponding top eigenvalue is 

T^((b) = v/7 + 2cos0 . 
The fourth-order truncation of the Hamiltonian, is 



^/3+2 cos 

^ 





-y/ 3+2 cos ( 

o 



3+8 cos 0+8 cos 2 < 
3+2 cos <f) 



3+8 cos 0+8 cos 2 < 
3+2 cos 



Its corresponding top eigenvalue is 



Tj 4 )(0) = V2, 



' 3 cos 2 



7 cos + 6 + a 



3 + 2 cos < 



where 



a = ^9 cos 4 + 26 cos 3 + 45 cos 2 + 54 cos + 27. 

In Fig. 2, we show the superconducting transition tem- 
peratures, 

AT< n »($) = T c (0) - T C W($), 

as functions of <E> = 0/2-7T for various values of n, for the 
square network obtained from the truncated Hamiltoni- 
ans H^ n \ Here T c (0) equals 4, which is the largest eigen- 
value of tight-binding electrons confined on the square 
lattice in the absence of a magnetic field. It is important 
to stress that as the order of approximation is increased, 
more geometrical information of the lattice is included 
in the interference treatment, and more fine structures 
are resolved. At every step, i.e., for a given size of the 
network, we can observe the corresponding dips appear- 
ing and then becoming sharper. We emphasize that our 
highest-order approximant, Tc 70 \&) has closely reached 
the infinite system size limit, AT C (<1>). The flux values 
where the cusps/dips occur have also been labeled. 



V. SELF-SIMILARITY IN THE PHASE 
BOUNDARY OF THE SUPERCONDUCTING 
SQUARE WIRE NETWORK 

In this section, we explicitly demonstrate an important 
property: the self-similarity of the phase boundary of the 
superconducting square wire network. This is exempli- 
fied in Fig. 3, where we use AT C (70) ($) for AT C ($) and 
omit the superscript. In (a), we plot AT C ($) for <!> in the 
interval between and 1. In (b) and (c), we plot AT c (<j>) 



for $, respectively, in the ranges [0.333 ~ 1/3, 0.4765] 
and [0.5235, 0.667 ~ 2/3]. Figures (b) and (c) can be re- 
garded as the first generation of the original diagram (a) , 
in the sense that (b) is enlarged from the maximum in 
the left part of (a) and (c) is enlarged from the maximum 
in the right part of (a). 

This enlargement process is continued as follows: (d) 
with $ e [0.375 = 3/8, 0.3978] and (e) with $ e 
[0.4025, 0.4286 ~ 3/7] are, respectively, the enlarge- 
ments of the left and right maxima of (b). Similarly, 
(f) with $ e [0.5714 ~ 4/7, 0.5975] and (g) with 
$ S [0.6022, 0.625 = 5/8] are, respectively, the enlarge- 
ments of the left and right maxima of (c). Figures (d), 

(e) , (f), and (g) can be regarded as the second genera- 
tion of the original phase diagram (a). In this way, it 
is straightforward to deduce that the third generation of 
(a) will consist of 8 phase diagrams: each of (d), (e), 

(f) , and (g) contributes two diagrams. It is evident that 
these phase diagrams resemble one another except that 
the phase diagrams gradually become asymmetric. 

As shown in these figures, we also label the values of 
$ indicating the cusps/dips in AT C (<1>). These nine flux 
values are characteristic of each phase diagram. Indeed, 
there are general relations between these sets of flux val- 
ues in different generations. Let {po/qo} represent the 
set of these flux values in (a), i.e., po/qo = 1/4, 2/7, 1/3, 
2/5, 1/2, 3/5, 2/3, 5/7, and 3/4. Denoting the set of 
characteristic flux values in any of the phase diagrams in 
the first generation by {pi/qi}, we find that the corre- 
sponding flux values in (b) are given by 



Pi 

qi 



qo 



3<?o - Po 
and those in (c) are given by 



Pi 



Pa + qo 



qi po + 2% 

For instance, given po/qo = 1/2 in (a), we have the cor- 
responding 

pi/qi = 2/(6-1) = 2/5 

in (b) and 

Pi/«i = (l + 2)/(l + 4) = 3/5 

in (c). Furthermore, let \jpilqi\ stand for the sets of 
the corresponding flux values in the second-generation 
diagrams. In the second-generation diagrams [(d)-(g)] 
only 5 characteristic cusps/dips out of 9 are observable. 
There we find that the 7J2/92 hi (d) are related to the 
pi/qi in (b) by 

P2 qi 

q 2 3qi - pi ' 

jj 2 /<?2 in (e) are related to pi/qi in (c) by P2/92 = 
?i/(3#i — P\), P2/92 in (f) are related to p\/q\ in (b) 

by 



El 

92 



Pi + qi 
Pi + 2gi ' 
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and pi/<l2 in (g) are related to pi/qi in (c) by p 2 /q2 = 
(pi + «i)/(pi+2«i). 

We now summarize our construction of the hierarchy 
of these phase diagrams. As discussed previously, every 
diagram can generate two diagrams to the next gener- 
ation: one is enlarged from the left maximum and the 
other from the right maximum of this diagram. Thus, 
starting from the original phase diagram, i.e., AT C (<I>) 
for $ £ [0,1], we can generate 2™ diagrams to the nth 
generation for n > 1. Furthermore, each diagram covers 
a distinct range of $ from $ m i n to $ max . Let us arrange 
these diagrams in the following way, as we did in Fig. 3. 
We put all the diagrams belonging to the same genera- 
tion in a row in such an order that from the left to the 
right $ m i n (or $ max ) increases from the smallest to the 
largest. It is evident that half of them (2"" 1 diagrams) 
have $ max < 1/2 and the other half have $ m i n > 1/2. It 
is not difficult to see that this kind of arrangement will 
be automatically satisfied in the following way. Following 
the same order of the diagrams in the previous generation 
and using them one by one, we put two new generated 
diagrams side by side with the one from the left maxi- 
mum to the left and the one from the right maximum to 
the right. It is interesting to notice that, for each genera- 
tion, the diagrams located at the left part of $ — 1/2 are 
mirror images of those located at the right part. This 
symmetry originates from the property that the phase 
diagram of AT C (4>) with 4> 6 [0, 1] is symmetric around 
$ = 1/2. 

Indeed, there are one-to-one correspondences between 
the sets of the characteristic flux values, where cusps/ dips 
in the phase boundaries occur, in different generations. 
Let us label the diagrams from left to right in the nth 
generation by with i running from 1 to 2™ . Similarly, 
the diagrams in the (n + l)th generation are labeled by 
V^ n+1} with i running from 1 to 2 n+1 . Now let {p n /q n } 
represent the sets of the flux values characterizing the 
cusps/dips in AT C ($) in any of the phase diagrams in 
the nth generation and {p n +\/q n +i\ be the sets belong- 
ing to the diagrams in the (n -I- l)th generation. The re- 
lations between the (p n+ i/<7„ + i)'s and the (p n /q n )'s are 
as follows. For 1 < i < 2™, the p n +i/?n+i in the diagram 
T)Y l+l > [one of the diagrams in the (n + l)th generation 
that located on the left hand side of $ = 1/2] is related 
to the p n /q n in X>,- by 

Pn+i _ q n 
q n +i 3q n -p n 

and for 2" + 1 < i < 2 n+1 , the p n +i/qn+i in the diagram 

v\ n+1) [the second half of the diagrams in the (n + l)th 
generation that located on the right hand side of <& = 1/2] 

is related to the p n /q n in T>^_ 2n by 

Pn+l _ Pn + q n 

q n +l Pn + 2<?„ 
Self-similarity in the AT C (<£>) curve is a consequence of 



the fractal energy spectrum of Bloch electrons in a mag=. 
netic field which was examined in detail by Hofstadter.EZl 
However, as far as we are aware, the explicit derivation 
of the self-similarity of the measurable part, the lowest 
energy state, was not presented before. 

Recently, the influence of classical chaos on Xhis so 
called "Hofstadter's butterfly" has been studiedEa Fur- 
thermore, a semiclassical theory for the dynamics of 
electrons in a magnetic Bloch band has been devel- 
oped and rused to explain the clustering structure of the 
spectrum.E3 

VI. HONEYCOMB LATTICE 

For the honeycomb lattice, we denote the lattice path 
integrals by h 2 i- In other words, h 2 i is the exact sum of 
the phase factors of all 2/-step closed paths on the hon- 
eycomb lattice. In this section, 0/2-7T corresponds to the 
magnetic flux through an elementary honeycomb plaque- 
tte, i.e., 

_ 3V3c 2 B 
2rr ~ 2$ 
The results for h-z, h±, . . . , /120 are: 

h 2 = 3, 
hi = 15, 

he = 87 + 6 cos 0, 

tig = 543 + 96 cos 0, 
h w = 3543 + 1080 cos + 30 cos 20, 
h 12 = 23859 + 10560 cos (j> + 726 cos 20 + 24 cos 30, 
h u = 164769 + 96096 cos 0+11130 cos 20 + 798 cos 30 
+42 cos 40, 

/lie = 1162719 + 839040 cos 0+ 138720 cos 20 
+ 15648 cos 30+ 1536 cos 40 + 96 cos50, 

/lis = 8363895 + 7143210 cos + 1537668 cos 20 
+237714 cos 30 + 33246 cos 40 + 3834 cos 50 
+252 cos 60+ 18 cos 70, 

h 20 = 61216275 + 59862000 cos0+ 15829200 cos 20 
+3103320 cos 30 + 555390 cos 40 + 89520 cos 50 
+ 10920 cos 60+ 1320 cos 70+ 120 cos 80. 

Notice that hi and /14 involve paths that enclose zero 
net flux. There are three closed paths of 2-steps each. 
Thus, /12 = 3, the coordination number of the lattice. 
/i6 is the first lattice path integral with a net flux (in 
this case flux through one hexagon). There are three 
counter-clockwise and three clockwise six-step paths go- 
ing through a hexagon. Thus, the term 6cos0 in he- 
It is possible to derive the first few path-integrals ana- 
lytically "by hand" , by just counting paths and keeping 
track of the enclosed flux. The longer-length ones can be 
computed via symbolic manipulation software. 
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We have computed the lattice path integrals for the 
honeycomb lattice up to /i206> which are obtained by ex- 
actly summing up ~ 10 96 closed paths. These calculated 
/i2i's have enabled us to obtain the phase boundary up 
toTi 1O4) (0). 

It is straightforward to compute the non-zero parame- 
ters b n from the obtained results for h^i . The correspond- 
ing truncated Hamiltonians, H^ n \ can then be readily 
constructed. For instance, the second-order truncation 
of the Hamiltonian is 



and 



H 



(2) 



V3 
V3 



Its corresponding top eigenvalue is = V3- The third- 
order truncation of the Hamiltonian is 



V3 
V3 V2 
V2 



Its corresponding top eigenvalue is = \/5. Both T c 

(3) 

and T c w; are independent of 0. This is understandable 
from the fact that the shortest length for a closed path 
on the honeycomb lattice to enclose the magnetic flux 
is for / = 6 while and only contain elements 
derived from /X2 and \x^. The fourth-order truncation of 
the Hamiltonian is 



(2) 



o Vs o o 
o o 

V2 V2" + COS ( 

V2 + cos 



Its corresponding top eigenvalue is 



In Fig. 4, we show the superconducting transition tem- 



peratures, Ar c (n) ($) = T c (0) - T c w ($), as functions of 
$ = 4>/2tt for various n for the honeycomb network ob- 
tained from the truncated Hamiltonians H^ n ' . Here T c (0) 
equals 3, which is the largest eigenvalue of tight-binding 
electrons confined on the honeycomb lattice in the ab- 
sence of a magnetic field. 

We observe that as the order of approximation is in- 
creased, more geometrical information of the lattice is in- 
cluded in the interference treatment, and more fine struc- 
tures are resolved. This explains the origin of the fine 
structure observed: the more geometric information on 
the lattice is explored by the paths of the electrons, the 
sharper the fine structures. 

We emphasize that our highest-order approximant, 
ri 104 '($) has closely reached the infinite system size 
limit, AT C ($). The flux values where the cusps/dips oc- 
curred have also been labeled. In general, besides the 
cusp at $ = 1/2, there are cusps at 

m 

$ = 



■,(n). 



m+ 1 
2m + 1 



with m > 1. Our computed phase boundary comi 
well with the observed cusps present in experimentsEJ' 21 



VII. TRIANGULAR LATTICE 

For the triangular lattice, we denote the lattice path 
integrals by t\. In other words, U is the exact sum of the 
phase factors of all l-step closed paths on the triangular 
lattice. In this section, 0/27T corresponds to the magnetic 
flux through an elementary triangular plaquette, i.e., 



2tt 



V3c 2 B 
4$ 



The results for t% through iio are: 



6, 

12 cos0, 

66 + 24 cos 2^, 

300 cos + 60 cos 30, 

1020 + 840 cos 20+168 cos 40+12 cos 60, 
6888 cos + 2604 cos 30 + 504 cos 50 + 84 cos 70, 
19890 + 23904 cos 20 + 8568 cos 40 + 1968 cos 60 
+432 cos 80 + 48 cos 100, 
164124 cos + 85944 cos 30 + 29628 cos 50 
+8496 cos 70+ 1980 cos 90 + 432 cos 110 
+36 cos 130, 

449976 + 654840 cos 20 + 317940 cos 40 

+ 1 14360 cos 60 + 37560 cos 80 + 10380 cos 100 

+2700 cos 120 + 540 cos 140 + 60 cos 160. 



Here we explain how the first few lattice path integrals 
are obtained. Since there is no path of one step for re- 
turning an electron to its initial site, t\ is always equal 
to zero. There are six closed paths of two steps each [re- 
tracing each other on one bond (■<-►), where the dot • 
indicates the initial site], thus 



t 2 

h 
h 
h 
t 6 

*7 



tg 



h = 6 



= 6 e l00 = 6 = 



where z is the coordination number of the lattice. 

There arc 12 three-step closed-paths enclosing a trian- 
gular cell [three counterclockwise (-v)i an d three clock- 
wise (*v)]- Thus 



6-V +6-V 



6e 



6 e 



12 cosi 



2m + 1 



There are 66 closed paths of four steps each enclosing 
zero flux each: six retracing twice on one bond (■«-►); 
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thirty starting from a site connecting two adjacent bonds 
and retracing once on each bond (<-> • W); and thirty 
moving two bonds away and then two bonds back to the 
original site (• — ► — ►). Since all of them enclose no area 
(i.e., no flux), then 



jnoflux = g 



30- 



30 



= 66. 



Among the 4-step closed paths, 24 of them enclose ad- 
jacent cells enclosing two triangles (12 counterclockwise 
and 12 clockwise) and contribute 



t 



two cells only 



12e 



2i<p 



12e 



-2i<j> 



24 cos 20 



to t^. Thus, it follows that t 4 — 66 + 24 cos 20. 

Note that t 2 i {t^i+i) consist of only even (odd) har- 
monics of the flux. We have computed the lattice path 
integrals for the triangular lattice up to tug, which are 
obtained by exactly summing up ~ 10 90 closed paths. 
These calculated ti 's have enabled us to obtain the phase 
boundary up to ri 6t ^(0). 

By using the calculated results for ti, the parameters 
a n and b n , and subsequently the corresponding truncated 
Hamiltonians, H^ n \ can be obtained. For instance, the 
second-order truncation of the Hamiltonian is 



H 



(2) 



V6 
V6 2cos<; 



Its corresponding top eigenvalue is 

Ti 2 )(0) 



cos + v6 + cos 2 0. 
The third-order truncation of the Hamiltonian is 
V6 



ff (3) 





V6 



2 cos< 



y/i +4 COS 2 



-8cos< 



■ 4 cos 2 
+ 16 cos 3 . 



1 + 4 cos 2 i 

-.(3), 



Its corresponding top eigenvalue Tc (0) can also be ob- 
tained analytically. 

In Fig. 5, we show the superconducting transition tem- 
peratures, AT c (n) ($) = T c (0) - T c (n) ($), as functions of 
<f> = (p/2n for various n for the triangular network ob- 
tained from the truncated Hamiltonians . Here T c (0) 
equals 6, which is the largest eigenvalue of tight-binding 
electrons confined on the triangular lattice in the absence 
of a magnetic field. The following physical picture is 
clear from those plots: as the order of approximation 
is increased, more geometrical information of the lattice 
is included in the interference treatment, and more fine 
structures are resolved. 

Our highest-order approximant, Tc 60 \&) has closely 
reached the infinite system size limit, AX" C (<I>). The flux 
values where the cusps occurred have also been labeled. 
In general, besides the cusps at $ — 1/2, 1/5, 4/5, 5/16, 
and 11/16, there are cusps/dips at 



$ = 



and 



with m > 1. 



<f> = 



2m + 2 



VIII. KAGOME LATTICE 

Our computed phase boundary compares well with 
the observeiL .cusps present in a series of interesting 
experimentfCJ'EJ. ■— ■ r— ■ ■— ■ 

For the kagome latticeJiSEJc-l we denote the lattice 
path integrals by fc;. In other words, ki is the exact 
sum of the phase factors of all /-step closed paths on the 
kagome lattice. Here (f>/2ir corresponds to the magnetic 
flux through an elementary triangular plaquette, i.e., 

j}_ _ V3c 2 B 
2?r _ 4$ 
The results for fc 2 through kn are: 

&3 = 4 COS (f), 

k 4 = 28, 

k 5 = 60 cos 4>, 

k 6 = 244 + 16 cos 2^ + 4 cos 60, 

k 7 = 756 cos + 28 cos 70, 

k 8 = 2412 + 416 cos 20 + 96 cos 60 + 80 cos 80, 

kg = 9216 cos + 76 cos 30 + 36 cos 50 

+756 cos 70+ 120 cos 90, 
fcio = 25804 + 7560 cos 20 + 1860 cos 60 

+2480 cos 80 + 100 cos 100 + 20 cos 140, 
fcn = 112420 cos + 2816 cos 30 + 1276 cos 50 

+ 14608 cos 70 + 4400 cos 90 + 44 cos 1 1 

+44 cos 130 + 176 cos 150. 



Note that k^i (fc2/+i) comprise only even (odd) harmonics 
of the flux. We have computed the lattice path integrals 
for the kagome lattice up to £gg, which are obtained by 
exactly summing up ~ 10 58 closed paths. These calcu- 
lated fc;'s have enabled us to obtain the phase boundary 
up toT c (5O) (0). 

By using the calculated results for ki, the parameters 
a n and b n , and subsequently the corresponding truncated 
Hamiltonians, H^ n \ can be obtained. For instance, the 
second-order truncation of the Hamiltonian is 



2 
2 cose 



2m + 2 



Its corresponding top eigenvalue is 

Tc 2) i<t>) = o ( cos + V16 + COS 2 0) . 
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The third-order truncation of the Hamiltonian is 

" 2 

2 cos 6 



y/3 - COS 2 



cos^ 



COS (j) + COS 

3 — cos 2 <b 



(3) 

Its corresponding top eigenvalue T c (cj>) can also be ob- 
tained analytically. 

In Fig. 6, we show the superconducting transition tem- 
peratures, AT c (n) ($) = T c (0) - T c (n) ($), as functions of 
<f> = 0/27T for various n for the kagome network obtained 
from the truncated Hamiltonians . Here T c (0) equals 
4 which is the largest eigenvalue of tight-binding elec- 
trons confined on the kagome lattice in the absence of 
the magnetic field. It is seen that as the order of the ap- 
proximation is increased, more geometrical information 
of the lattice is included in the interference treatment, 
and more fine structures are resolved. We emphasize 
that our highest-order approximant, Tc 50 \&) has closely 
reached the infinite system size limit, AT c (<j>). The flux 
values where the cusps/dips occurred have also been la- 
beled. Our computed phase boundary compares well 
with the obserjied cusps present in a series of interesting 



expe 
Rcf.t 



imentscjO. See also the systematic calculations in 



IX. DISCUSSION 

In the following, we discuss the general trends in the 
approximants for these phase diagrams presented in the 
above sections. 



A. Comparison of the structure in the phase 
boundaries 

In the lower-order approximants, the first noticeable 
development in the phase boundaries of square, honey- 
comb, and triangular lattice is the formation of dips when 
the flux per elementary plaquette is equal to to$o/2, 
where m is an integer. When the order of approximation 
is increased, the dips at $ = 1/2 become sharper and at 
the same time more fine structures (other local minima) 
begin to emerge. Eventually, the dips at various different 
flux values become cusps. 

It is interesting to notice that, among these three lat- 
tices, the development of the cusps is most rapid for the 
triangular case while the honeycomb is slowest. This dif- 
ference originates from the fact that for identical length, 
lattice path integrals for the triangular lattice contain the 
richest quantum interference effects because the number 
of paths and the areas they enclose are both the largest. 
For the kagome network, the rapid development of cusps 
at 

$ kagomo = 1 I 3 5 3 7 
8' 4' 8' 8' 4' a M 8 



can be seen from lower-order approximants. For addi- 
tional discussion on the kagome case, see RcfM For ex- 
tensions of these techniques to other problems, see Ref.Ea. 

In general, the resulting phase diagrams — with the oc- 
currence of cusps/dips at different sets of flux values — 
are a direct consequence of the geometries of the lat- 
tices, which is explicitly reflected in the corresponding 
expressions of the lattice path integrals. We stress that 
our evaluation of the lattice path integrals to extremely 
long lengths has enabled our calculated T C (B) to achieve 
close convergence to the infinite system size. Indeed, for 
n ~ 10, important features in the phase boundaries of 
square, triangular, and kagome networks are well devel- 
oped. 

Finally, in order to facilitate a comparison between the 
different phase boundaries, in Fig. 7 we plot AT C (<!>) as a 
function of for the square, honeycomb, triangular, and 
kagome superconducting networks. Here the AT c (<l>)'s 
are taken from their respective highest-order approxi- 
mants and 4> is the flux through their respective elemen- 
tary cells as discussed in the previous sections. Here we 
omit the subscripts indicating the order of approxima- 
tion. The values of the magnetic flux corresponding to a 
number of prominent cusps/dips are also labeled. 



B. Comparison of the phase boundaries of the 
single-loop and lattice cases 

From Figs. 1 and 7 (a-c), we can readily see the dif- 
ferences between the phase boundaries of a single super- 
conducting cell and its corresponding superconducting 
network. For both cases, AT C vary periodically with 
the magnetic flux through a single elementary cell and 
have the same period $o of oscillation. We now focus on 
AT C ($) for $ in the interval between and 1. AT C ($) 
is symmetric around $ = 1/2. However, there are many 
distinct features between AT C of a single cell and that of a 
network. These differences are due to long-range correla- 
tions of the many-loops effect present in the lattice. For a 
single superconducting cell, AT' c (<i>) increases monotoni- 
cally from $ = to $ = 1/2 and decreases monotonically 
from $ = 1/2 to $ = 1. The maximum at $ = 1/2 ex- 
hibits a sharp peak. Indeed, the overall shape of AT C (&) 
resembles the combination of two identical half parabolas 
both reaching their maximum at (f> = 1/2. On the con- 
trary, the overall shape of AT C (Q) for the corresponding 
superconducting networks look like downward parabolas 
with many local cusps between $ = and $ = 1. The 
most prominent cusps are located at $ = 1/2. The posi- 
tions of other cusps/dips depend on the underlying lattice 
types of the networks. 
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C. Differences between our approach and the 
traditional moments and Lanczos methods 

In electronic structure calculations there is a method 
to compute the density of states called the moments 
method. This is similar to our approach in the sense 
that hi can be interpreted as the moments, (\E'j|/f'|\E , j). 
However there are several important differences between 
the standard "moments method" and our problem. The 
typical use of the moments method: (i) focuses on the 
computation of electronic density of states (instead of 
superconducting T c 's); (ii) is totally numerical (instead 
of mostly analytical); (iii) is done at zero magnetic field 
(instead of obtaining expressions with an explicit field 
dependence); (iv) does not focus on the explicit compu- 
tation of lattice path integrals; and (v) does not study 
the physical effects of quantum interference (which is at 
the heart of our calculation and physical interpretation) . 
In conclusion, the traditional use of the moments method 
in solid state is significantly different from the approach 
and problem studied here. 

Another way to diagonalize Hamiltonians is called the 
Lanczos method. This method directly obtains the tridi- 
agonal form, without computing the moments; thus dif- 
fers in a significant way from the approach used here 
(where the explicit computation of the moments is one 
of our goals, since they can be used for other electronic 
property calculations). Furthermore, it is not convenient 
to use standard Lanczos method in our particular prob- 
lem because it is extremely difficult to directly derive the 
parameters and the states of the iterative tridiagonaliza- 
tion procedure. This is so because of the presence of the 
magnetic field. On the other hand, the moments method 
provides standard procedures to diagonalize a matrix af- 
ter the moments are computed. 

D. Commensurability and other Matching Effects 

An essential physics issue in this problem is commensu- 
rability. Another one is quantum interference — due to the 
motion of electrons in multiconnected geometries. This 
section briefly overviews related systems where commen- 
surability and matching effects (due to externally applied 
magnetic fields) play an important role. The first exam- 
ple will be flux pinning. 

Flux pinning in type II superconductors is of both 
technological and scientific interest. While most exper- 
iments focus on the effects of random pinning distribu- 
tions, some investigations have been carried out on peri- 
odic arrays of pinning siteso. These find striking peaks in 
the magnetizatior£3 and critical current J c . These peaks 
are believed to arise from the greatly enhanced pinning 
that occurs when parts of the vortex lattice become com- 
mensurate with (i.e., match) the underlying periodic ar- 
ray of pinning sites. Under such conditions, high-stability 



vortex configurations are produced which persist under 
an increasing current or external field. 

Other important vortex matching effects have also 
recently been obse:cued-| in a variety of different super- 
conducting systemsc3E9, including long-Josephson junc- 
tions with periodically-spaced groovesLJ, superconduct- 
ing networkscj, and the matching of the VL to thexrystal 
structure of YBa2Cu30y due to intrinsic pinningEH 

Matching effects between a vortex lattice and periodic 
pinning arrays produce a rich variety of effectsEil. The dy- 
namics observed in these systems is quite different from 
the one faund for random arrays of pinning sites (see, 
e.g., Ref.E3 and references therein). 

Non-superconducting systems also exhibit magnetic- 
field-tuned matching effects, notably in relation to elec- 
tron motion in periodic structures where unusual behav- 
iors arise due to the incommensurability of the magnetic 
length with the lattice spacing. A recent example of these 
is provided by the anomalous Hall plateaus of "electron 
pinbaH"E2l orbits scattering from a regular array of anti- 
dots. 

Commensurate effects also play central roles in many 
other areas of physics, including plasmas, nonlinear 
dynamics^, the growth of crystal surfaces, domain walls 
in incommensurate solids, quasicrystals, Wigner crystals, 
as well as spin and charge density waves. The next sec- 
tion discusses in some detail an example in nonlinear dy- 
namics (which is virtually unknown in the solid state lit- 
erature) that produces a fractal phase boundary which is 
strikingly similar to the one measured for square super- 
conducting networks — because both are determined by 
commensurability effects. 

E. Kagome-pinned vortices: "Correlated Melting" 
and Cooperative Ring Excitations for 
Doubly-Degenerate Ground States 

Notice that the fluxoid configurations for f — 1/2 for 
the superconducting networkSj-fe-g-, Fig. 3 of Yi Xiao, 
et al, in the companion articleo) has two ground states 
that correspond to the two degenerate ground states of 
the second matching field of vortices in type II supercon- 
ductors with a kagome periodic array of pinning—sites. 
The latter has been systematically studied in Ref.Ej. 

The kagome pinning potential at the second match- 
ing field shows novel .and interesting dynamics as a func- 
tion of temperature,Ej including a phase with rotating 
vortex triangles caged by kagome hexagons ("coopera- 
tive ring elementary excitations"), and there is geometric 
frustration for T —> with a doubly-degenerate ground 
state. At finite temperatures, the three vortices inside 
the kagome' hexagon can move and rotate by 60 degrees. 
This is done cooperatively by the three vortices. They 
motion is similar to the "cooperative ring exchange" mo- 
tion proposed by Feynman for elementary excitations in 
Helium 4. 
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In other words, for the second matching field for the 
kagome pinning lattice, the elementary excitation of the 
three interstitial vortices is a 60 degree rotation, rotating 
as a cooperative ring. These type of collective or corre- 
lated cooperative ring exchange has also been studied in 
the context of the quantum Hall effect. r— ■ 

For increasing temperatures, a novel type of meltingo 
appears. This can be described as "correlated melting" in 
the sense that the "triangle" or "loop" first melts in the 
angular coordinate, while the radial coordinate does not 
melt until much higher temperatures are reached. The 
elementary excitations are the thermal analog of certain 
types of squeezed states (where fluctuations strongly af- 
fect a coordinate and less the other coordinate). They 
are also analogs of the rotational isomers or "comfor- 
mations" that are often found in molecules, where three 
atoms and molecules can cooperatively oscillate back and 
forth between two degenerate ground states. 

This tope of "controlled melting" or "correlated 
melting" Ea of the particles inside a potential energy trap 
could also be visualized with a colloidal suspension sur- 
rounded by six pinned (e.g., by laser tweezers) charged 
particles. This type of "vortex-analog" experiment is eas- 
ier to visualize (e.g., via optical microscope) r-than using 
vortices. Still, Lorentz microscopy techniques^ could di- 
rectly image such motions in the vortex case. 

F. Fractal Phase boundaries and Fractal Boundaries 
of Basins of Attraction 

There is a striking similarity between two apparently 
unrelated problems: the superconducting-normal phase 
boundary of a square superconducting network (our 
Fig. 2), and the fcactal phase boundary (see, for instance, 
Fig. 6.26 of Ref.Efl) of basins of attraction of a dynami- 
cal systems map studied last century by Weierstrass and 
generalized much later by Hardy in 1916. 

The reason for this very interesting similarity among 
these two apparently unrelated problems is because the 
commensurability condition dominates both problems 
and produces a large dip at 1/2, and smaller dips at 1/4, 
1/3, 2/5, etc., as discussed previously in this work. 

It is interesting to summarize how to obtain the Weieu- 
strass fractal boundary of two basins of attractions. 
Consider the dynamical map M 

(xfc+i, 6 k +i) = M(x k , k ) 

defined by, 

Xk+i = Ax fe + cosflfe 

and 

6 k+1 = 26 k (mod 2?r) . 

When 1 < A < 2, the map M has two attractors, at 
x = ±oo. Indeed, since the eigenvalues of the Jacobian 



matrix are 2 and A > 1, there are no finite attractors. 
Therefore, 

M N {x , 6 ) = (x Nl 6 N mod(27r)) , 

and x 7v tends to either +oo or — oo as N — + oo, except 
for the unstable boundary set 

* = f(0) , 

for which xn remains finite. 

To locate this x = f(0) boundary set, first note that 

6 k = 2 k 6 (mod 2tt) . 

The map is non-invertible since it is two-to-one. However, 
any xn can be selected and then find one orbit that ends 
at (xn, &n), by using the above k and taking 

Xk-i = A _1 [x fe - cos(2 k ~ 1 6 )] • 

For a given (xn, &n), this orbit starts at 

N-l 

x = \~ 1 x N - ^2 X ~ l ~ 1 cos(2'6» ) . 

1=0 

Those (atoj Oo) such that xn is finite as N — > oo, define 
the boundary x — f(6) between the two basins. There- 
fore the relation between these x and 6 is given by 

oo 

x = -Y / ^ 1 cos(2 l 6 ) = f(6) . (8) 

1=0 

This sum obviusly converges, since A > 1. However, 
its derivative 

1=0 v 7 

diverges, since A < 2. Thus, f(6) is non-differentiable. 
Like our superconducting-normal phase boundary, it has 
a large cusp at 1/2, and smaller cusps at 1/3, 1/4, 2/5, 
etc. Moreover, it is also symmetric around 1/2, and it 
strongly resembles the AT C ($) (obtained near R(T) = 0) 
for a square lattice. Indeed, AT c (<i>) corresponds to x, 
and $ to 0. 

The fractal dimension of x — f(0), Eq.(8) above, is 

, „ lnA 

dc = 2 - — . 
In 2 

The precise value of d c depends on the value of A. Re- 
call that 1 < A < 2. For A sligthly less than two, the 
fractal dimension d c approaches one, and the dips are 
not pronounced. This is similar to the superconducting- 
normal phase boundaries measured not too close to T c 
(e.g., at mid-point drop for the R(T) plot). When the 
phase-boundary is measured very near T c (when R(T) 
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is very near zero) , the number of discernible dips grows 
and they become very sharp (see, e.g., Figs. 10 and 11 
of Ref Ji-3) . This would corespond to A slightly above one; 
thus, the fractal dimension d c of the Weierstrass func- 
tion would be closer to two (i.e., a "rougher" or "spikier" 
curve). |— - 

Indeed, RefM solved for AT C ($) beyond the mean field 
theory approximation, obtaining a phase boundary sim- 
ilar to the Weierstrass function for A slightly above one, 
and d c near two--. That superconducting-normal phase 
boundary in ref.L3 has very sharp cusps and dips, and 
(like the Weierstrass function) it is a phase boundary 
between attractors. The map for the superconducting 
networks is obtained from a real-space renormalization- 
group technique. The mean field limit provides smoother 
phase boundary with A closer_tp one. The real-space 
bond-decimation scheme of Ref.liJ also favors fluxes of the 
form $ = 2 l $o- This is clear from the way the real-space 
renormalization-group scheme is constructed, where four 
elementary cells are "bio eked- away" into a larger cell 
with new renormalized effective couplings. Four of these 
super-cells are then blocked away into another, larger 
cell, enclosing 16 elementary cells (or 4 supercells). This 
process is iterated, until the renormalization group pro- 
cedure coverges (at the phase boundary) or diverges to 
fixed points located away from the fixecL-noint (e.g., +oo). 
This (beyond-mean-field) RG iterations and the Weier- 
strass iteration involve very similar types of maps and 
this generates the strikingly similar curves. 

In summary, the Weierstrass function and our real- 
space renormalization group approachtil both produce 
phase boundaries which are strikingly similar. In partic- 
ular, both are non-differetiable, symmetric around 1/2 
and have a very similar hierarchy of cusps. 

X. COMPARISON OF THE PHASE 
BOUNDARIES OF SUPERCONDUCTING 
HONEYCOMB AND KAGOME NETWORKS 

Here we discuss an interesting relation between the 
phase boundaries of superconducting honeycomb and 
kagome networks which is due to the geometrical arrange- 
ments of these two types of lattices. Indeed, and as kindly 
pointed out to us by Y. Xiao and P.M. Chaikin, it is very 
useful to focus on the region < $ < 1 for the honey- 
comb network and the region < $ < 1/8 for the kagome 
network. 

As shown in Fig. 8(a) through 8(d), though the over- 
all shapes of the phase diagrams are different, there is a 
one-to-one correspondence between the dips in the hon- 
eycomb AT C (<I>) for $, the flux through an elementary 
hexagon, in the range [0, 1] and those in the kagome 
AT C (<I>) for $, the flux through an elementary triangle, 
in the ranges [0, 1/8], [1/8, 1/4], and [1/4, 3/8]. To 
state this relationship more precisely, let {p/q} be the 
set of flux values characterizing a number of dips in the 



AT C ($) curve for the honeycomb network. For instance, 
as labeled in (a), 

{p/q} = 1/3, 2/5, 3/7, 1/2, 4/7, 3/5, and 2/3. 

It is observed that the corresponding set of flux values 
for the dips to occur in the kagome AT C ($) curve would 
be {p/8q} when $ lies in the range [0, 1/8]. Similarly, 
the corresponding sets read, respectively, 

r 1 P P + 9i 
* = {- + — = - — -} 

X 8 8q 8q 5 

for $ € [1/8, 1/4] and 

1 ^ = p + 2q 
4 8q 8q 5 

for $ e [1/4, 3/8]. Note that for $ in the range 
[1/4, 3/8], the dips in the AT C (<I>) curve become less evi- 
dent: only five flux values are observed and labeled. The 
location and magnitude of the dips found here are con- 
sistent with recent very interesting experiments by the 
NEC and Princeton groupsEJEII. 

Recall that kagome magnets are known to have de- 
generate ground states (see, e.g., Ref.E3 and references 
therein). Likewise, for superconducting kagome networks 
at half filling, there are several possible ways to arrange 
fiuxeSynproducing a large degeneracy in the T = ground 
state&3. This issue of degeneracy between two states has 
been systematically studied as a function of temperature 
via computer simulations on superconducting samples 
with a kagome-arranged periodic array of pinning sitescj. 
The second matching field in this system has two fluxons 
per pinning site. This corresponds to the / = 1/2 state in 
the kagome superconducting network. For this value of 
the externally applied magnetic field, every hexagon has 
two states (with entropy fee log 2). N hexagons would 
have 2 N states, and a very large entropy 

S (N hexagons) _ N k B \ og2 . 

Thus, at the second matching field, superconductors 
with either a kagome or an hexagonal array of pinning 
sites both have "low-energy states" with a very large 
degeneracy and a huge (low-T) entropy. Thus, when 
cooling from high temperatures, it is difficult to find a 
unique T = ground state. Transport measurements 
and mean field theory perhaps might not be sufficient 
to fully elucidate the role of bistability and degeneracy 
in this system. In order to explore this scenario in a 
more systematic, manner, different tools (e.gp-.flux imag- 
ing techniques^ and computer simulationscJ of vortex 
dynamics on kagome lattices) might be needed. 

After this work was completed, we became aware of 
a vexy interesting relevant work by Park and Huse in 
Ref.E3. Using Ginzburg-Landau theory, they study su- 
perconducting kagome wire networks in a transverse 
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magnetic field when the magnetic flux through an ele- 
mentary triangle is a half of a flux quantum. They cal- 
culate the helicity moduli of each phase to estimate the 
Kosterlitz-Thouless (KT) transition temperatures. At 
the KT temperatures, they estimate the barriers to move 
vortices and the effects that lift the large degeneracy in 
the possible flux patterns. 
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XI. SUMMARY 

In conclusion, we present a detailed study of the mean- 
field superconducting-normal phase boundaries of super- 
conducting square, honeycomb, triangular, and kagome 
networks. Our investigations are based on studying the 
quantum interference effects arising from the summation 
of all the closed paths the electron can take on the under- 
lying lattices. Other problems^ have also been studied 
in terms of quantum interference of electron paths. We 
then adopt a systematic approximation scheme, to obtain 
the spectral edges of the corresponding eigenvalue prob- 
lems, and relate the features in the phase boundaries with 
the geometry of the underlying lattice being explored by 
the moving electrons. When the electrons are allowed 
to explore a sizable region of the network, our calcula- 
tions have quickly reached very close convergence to the 
infinite system size results. There are two particular ad- 
vantageous aspects of this approach. First, it enables us 
to evaluate the superconducting transition temperature 
as a continuous function of the applied magnetic field. 
Second, it enables us to achieve a step-by-step derivation 
of the physical origin of the many structures in the phase 
diagrams — in terms of the regions of the lattice explored 
by the electrons. In particular, the larger the region of 
the network the electrons can explore (and thus more 
paths are available for the electron), the finer structure 
appears in the phase boundary, and the sharper the cusps 
become. We find many new interesting features in these 
phase diagrams, which compare well with experiments. 
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FIG. 1. The oscillatory phase boundary, AT C ($), for a sin- 
gle superconducting loop. The top curve corresponds to a 
triangle (dashed line), the middle a square (dotted line), and 
the bottom a hexagon (solid line). $ is the magnetic flux 
through these cells in units of $o. 

FIG. 2. Superconducting transition temperature for the 
square network as a continuous function of the applied mag- 
netic field: AT c (n) ($) = T c (0) - T c (n) ($) versus $, the mag- 
netic flux through an elementary square cell. In (a) we show 
the superconducting-normal phase boundaries computed from 
the truncated Hamiltonians, H^ n \ for <£> in the range between 
0.2 and 0.8. We omit the parts of AT c (n) ($) for $ G [0, 0.2] 
and [0.8, 1] since there are no interesting features in these 
portions of ATj" (<&). From top to bottom, the orders of 
truncation are n = 5 (top curve), 6, 7, 8, 10, 15, 23, 39, 
and 70. Note the development of fine structures and cusps. 
The convergence is monotonia Note also that the close- 
ness between the curves for AT C (39) ($) and AT C (70) ($) im- 
plies that AT C ' 70 ' ($) has achieved close convergence to the 
infinite system size AT C ($). The inset schematically depicts 
a square lattice. In (b), we plot AT C ($) for $ G [0.2, 0.8] 
and label the values of the magnetic flux where observable 
cusps/dips occur. They include $ = 1/4, 2/7, 3/10, 1/3, 
3/8, 2/5, 3/7, 1/2, 4/7, 3/5, 5/8, 2/3, 7/10, 5/7, and 3/4. 
Here AT C ($) = AT C (70) ($) = T c (0) - T c (70) ($), our calculated 
highest-order approximant. 



FIG. 3. Field-dependent transition temperature, AT C (<&), 
of the superconducting square network for various dif- 
ferent ranges of $: from (a) to (g), respectively, 
$ G [0,1], [0.333 ~ 1/3, 0.4765], [0.5235, 0.667 ~ 2/3], 
[0.375 = 3/8, 0.3978], [0.4025, 0.4286 ~ 3/7], 
[0.5714 ~ 4/7, 0.5975], and [0.6022, 0.625 = 5/8]. It is clear 
that (b) is enlarged from the maximum in the left part of (a) 
and (c) is enlarged from the maximum in the right part of 
(a). Similarly, (d) and (e) are, respectively, the enlargements 
of the left and right maxima of (b) while (f) and (g) are, re- 
spectively, the enlargements of the left and right maxima of 
(c). We also include the labeling of the values of $ where 
there are cusps/dips in AT C ($). For the relations between 
these sets of flux values in different frames, see the text. The 
self-similarity in the phase boundary can be concluded from 
the resemblance of these figures though an asymetry in the 
height develops in each successive magnification. 

FIG. 4. Superconducting transition temperature for the 
honeycomb network as a continuous function of the applied 
magnetic field: AT c (n) ($) = T c (0) - T c (n) ($) versus $, the 
magnetic flux through an elementary hexagonal cell. In (a) 
we show the superconducting-normal phase boundaries com- 
puted from the truncated Hamiltonians, _ff' n ', for $ in the 
range between 0.3 and 0.7. We omit the parts of AT c ' n '($) 
for <£> G [0, 0.3] and [0.7, 1] since there are no interesting 
features in these portions of ATj™'($). From top to bot- 
tom, the orders of truncation are n = 9 (top curve), 10, 13, 
16, 21, 31, 41, and 104. Note the development of fine struc- 
tures and cusps. The convergence is monotonic. We believe 
that ATi 104 ^($) has achieved close convergence to the infi- 
nite system size AT C ($). The inset schematically depicts a 
honeycomb lattice. In (b), we plot AT C (<I>) for $ G [0.3, 0.7] 
and label the values of the magnetic flux where observable 
cusps/dips occur. They include $ = 1/3, 2/5, 3/7, 4/9, 5/11, 
6/13, 7/15, 8/17, 1/2, 9/17, 8/15, 7/13, 6/11, 5/9, 4/7, 3/5, 
and 2/3. Here AT C ($) = AT C (104) ($) = T c (0) - T c (104) ($), 
our calculated highest-order approximant. 
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FIG. 5. Superconducting transition temperature for the 
triangular network as a continuous function of the applied 
magnetic field: AT c (n) ($) = T c (0) - T c (n) ($) versus $, the 
magnetic flux through an elementary triangular cell. In (a) 
we show the superconducting-normal phase boundaries com- 
puted from the truncated Hamiltonians, H^ n \ for $ in the 
range between 0.15 and 0.85. We omit the parts of AT c (n) ($) 
for $ 6 [0, 0.15] and [0.85, 1] since there are no interesting 
features in these portions of ATc™\&). From top to bot- 
tom, the orders of truncation are n = 5 (top curve), 6, 7, 
10, 15, 29, and 60. Note the development of fine structures 
and cusps. The convergence is monotonic and rapid. Note 
also that the closeness between the curves for AT C (29) ($) and 
AT C (60) ($) implies that AT C (60) ($) has achieved close conver- 
gence to the infinite system size AT C (<I>). The inset schemat- 
ically depicts a triangular lattice. In (b), we plot AT C ($) for 
"I? 6 [0.15, 0.85], our calculated highest-order approximation 
to AT C ($), and label the values of the magnetic flux where 
observable cusps/dips occur. They include $ = 1/5, 1/4, 
5/16, 1/3, 3/8, 2/5, 5/12, 3/7, 7/16, 4/9, 9/20, 1/2, 11/20, 
5/9, 9/16, 4/7, 7/12, 3/5, 5/8, 2/3, 11/16, 3/4, and 4/5. 
Here AT C (<I>) = AT C (60) (<I>) = T c (0) - T c (60) ($), our calculated 
highest-order approximant. 

FIG. 6. Superconducting transition temperature for the 
kagome network as a function of the applied magnetic field: 
AT c (n) ($) = T c (0) - T c (n) ($) versus $, the magnetic flux 
through an elementary triangular cell. In (a) we show the su- 
perconducting-normal phase boundaries computed from the 
truncated Hamiltonians, H^ n \ for $ in the range between 
and 1. ^From top to bottom, the orders of truncation are 
n = 4 (top curve), 5, 6, 8, 10, 19, and 50. Note the de- 
velopment of fine structures and cusps. The convergence is 
monotonic. Note also that the closeness between the curves 
for AT C (19) ($) and AtJ 5()) ($) implies that AT C (50) ($) has 
achieved close convergence to the infinite system size AT C (<E>). 
The inset schematically depicts a kagome lattice. In (b), wc 
plot AT C (<E>) for <E> G [0, 1] and label the values of the mag- 
netic flux where observable cusps/dips occur. They include 
$ = 1/12, 1/8, 4/25, 1/4, 1/3, 3/8, 5/8, 2/3, 3/4, 19/24, 7/8 
and 11/12. Here AT C ($) = AT C (50) ($) = T c (0) - T c (50) ($), 
our calculated highest-order approximant. Note the absence 
of the cusp at $ = 1/2. This distinct feature is in sharp 
contrast to the cases for square, honeycomb, and triangular 
networks. 

FIG. 7. AT c ($)'s as functions of 3> between and 1 for the 
superconducting square, honeycomb, triangular, and kagome 
networks, respectively, from (a) to (d). Notice the difference 
in the vertical scales. 

FIG. 8. AT C ($) versus <J?. (a) is for the superconducting 
honeycomb network for 3> in the range [0, 1]. (b), (c), and (d) 
are for the superconducting kagome network for respec- 
tively, in the ranges [0, 1/8], [1/8, 1/4], and [1/4, 3/8]. 
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